Charge transport in nanochannels: a molecular theory 
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We introduce a theoretical and numerical method to investigate the flow of charged fluid mixtures 
under extreme confinement. We model the electrolyte solution as a ternary mixture, comprising two 
ionic species of opposite charge and a third uncharged component. The microscopic approach is 
based on kinetic theory and is fully self-consistent. It allows to determine configurational prop- 
erties, such as layering near the confining walls, and the flow properties. We show that, under 
appropriate assumptions, the approach reproduces the phenomenological equations used to describe 
electrokinetic phenomena, without requiring the introduction of constitutive equations to determine 
the fluxes. Moreover, we model channels of arbitrary shape and nanometric roughness, features that 
have important repercussions on the transport properties of these systems. 

Numerical simulations are obtained by solving the evolution dynamics of the one-particle phase- 
space distributions of each species by means of a Lattice Boltzmann method for flows in straight 
and wedged channels. Results are presented for the microscopic density, the velocity profiles and for 
the volumetric and charge flow-rates. Strong departures from electroneutrality are shown to appear 
at molecular level. 

I. INTRODUCTION 

Understanding the motion of electrical charges in liquid solutions is at the core of many recent advances in bio- 
logical and technological applications, ranging from molecular separation processes, DNA sequencing, oil recovery, 
decontamination of soils by extracting and removing pollutants and toxins, manipulation of colloidal materials, and 
so forth. Electro-nanofluidics allows controlling the movement of liquids in and around objects with characteristic size 
H of the order of 10 — 50 nm, and require a deep understanding of the transport properties of small amounts of fluids 
[1-3 . In fact, charged liquid solutions under confinement behave in a distinctive way under the influence of electric 
fields, because their flow properties depend on the charges of the moving ions and on those sitting at the pore walls. 
An ionic atmosphere oppositely charged to the surface charge and of thickness proportional to the Debye length, Xrj, 
forms next to the walls, giving rise to the so-called electric double layer (EDL). The difference with neutral fluids is 
more pronounced when the Debye screening length of the mobile charges becomes comparable with the size of the 
confining pores |4j [5] . 

Smoluchowski provided the first explanation of mass transport of electrolytes in electro-osmotic conduits by using 
the theory of the EDL of Gouy and Chapman. He computed the electric field due to the surface charges and to 
the mobile ions by solving the Navier-Stokes (NS) equation in the limit of small Reynolds numbers. In electrically 
driven flows, the predicted fluid velocity profile of an electrolyte displays a characteristic plug-like profile departing 
from the Poiseuille parabolic profile observed in pressure driven fluids [5|. The work of Nernst and Planck, instead, 
extended the standard theory of diffusion to include ionic systems, by representing the ionic current in terms of a 
phenomenological convection-diffusion-migration equation [7]. 

Although successful in describing microfluidic phenomena, the classical electrokinetic approach fails to represent the 
behavior of nanoscopic system. Its major shortcoming is the continuum description, which breaks down at atomistic 
scales, where the system behavior becomes dominated by wall effects and it is even hard to define a bulk volume. 
While in microfluidics the effect of the surface charge affects the fluid only within a distance of the order of \e> << H . 
in nanochannels the Debye length may cause the channel to become ion-selective by excluding one type of ion over the 
other [8j. Therefore, some features of the continuum approach have to be modified a) by considering the molecular 
nature of the solvent which competes for room in the pores and affects the structural properties of the EDL, and b) 
by treating the transport properties beyond the standard NS and Poisson-Nernst-Planck (PNP) equations which are 
valid in the hydrodynamic limit [5]. 

A variety of methods in statistical mechanics have been applied to study the transport and structural properties 
of charged fluids from a microscopic perspective. These methods include Molecular Dynamics (MD) simulations |10j . 
Density Functional theory in conjunction with the Navier-Stokes equation |llj . and an energy variational method 
|12| . Among these, the nonequilibrium MD technique [13 provides in principle the best description being subject to 
fewer approximations than other methods. However, in the case of electrokinetic flows it has to cope with two kinds 
of difficulties: a) in the systems of interest the ion concentrations can be orders of magnitude lower than the solvent 
concentration, so that statistical accuracy is reached only at the expenses of long computer runs, b) the electro-osmotic 
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speeds are much smaller than the thermal speed, rendering difficult to extract the signal from the background noise. 

The density functional approach of Griffith and Nilson deals rather accurately and microscopically with the 
configurational description of a ternary system describing an assembly of positive and negative ions immersed in 
a sea of hard sphere dielectric particles mimicking water molecules, but resorts to a macroscopic treatment of the 
fluid velocity introduced by means of a Navier-Stokes equation, thus neglecting the local dependence of the transport 
coefficient on the density fluctuations. The energy variational treatment of ions in protein channels is similar in spirit, 
although derived from different premises, but considers only a primitive model of electrolyte, that is, the solvent 
is not dealt explicitly |12j . For complex pore geometries of physical interest one usually resorts to numerical tools 
such as finite element method, finite difference, boundary element j 14] . spectral techniques |15| and recently Lattice 
Boltzmann methods [16-21], to mention just a few. 

In the present paper we propose to treat structural and dynamical aspects of electrokinetics in confined ionic systems 
by an extension of the non-local Boltzmann-Enskog equation |22| recently developed by us for neutral multicomponent 
mixtures. The method requires the knowledge of the inhomogeneous microscopic pair correlation functions of the fluid 
and in principle accommodates the same level of microscopic details as the Dynamical Density Functional theory [231 - 
125] . The advantage of the present approach is that both the equilibrium and the non equilibrium forces (which 
determine the transport coefficients) are derived from the same principle. Thus, we are able to predict not only the 
density profiles of the each species, but also the velocity profiles, without resorting to further approximations. In 
practice, the method brings together molecular aspects (microscopic structure, non-ideal gas effects in the equation 
of state, microscopically determined transport coefficients) with hydrodynamic aspects of correlated fluids. 

Our theory employs an approximate method, as introduced about twenty years ago by Dufty and coworkers |26il27|. 
in order to reduce the complexity of the full collision operator to a form amenable to fast iterative numerical solutions. 
This approximation is based on the idea of separating the fast degrees of freedom, giving a minor contribution to the 
transport properties, from the slow hydrodynamic modes, which on the contrary are very relevant. By this strategy one 
achieves a considerable simplification of the resulting kinetic equations with respect to the original Enskog-Boltzmann 
equations. This set of equations represents the interactions as self-consistent non-local forces and are amenable to 
numerical solution. In the past, numerical solutions of the Boltzmann-Enskog equation for inhomogeneous systems 
were prevented by the difficulty of solving for the distribution function in the 7-dimensional phase-space. This 
obstacle, however, can be overcome by following the same strategy introduced by us to handle neutral systems under 
inhomogeneous conditions. We have shown that a numerical solution method can be obtained in the framework of 
the Lattice Boltzmann approach [28-32 . 

The structure of this paper is the following: in section |TT] we prepare the scene, by briefly summarizing the phe- 
nomenological transport equations to the mass and charge transport in pores. In section |TTT| we set up the necessary 
formalism, based on an extension of the Enskog-Boltzmann transport equation, to deal with a ternary mixture com- 
prising positively and negatively charged hard spheres to simulate the the ions and uncharged hard spheres for the 
solvent. Next, by projecting the transport equations onto the hydrodynamic space we formally recover the same 
structure of equations for the charge, the mass and the momentum as the phenomenological model. We then obtain 
a microscopic expression of all the transport coefficients and include non ideal gas effects into the theory. We extract 
qualitative predictions from our theory and finally present solutions of the full coupled system of transport equations 
using the LBE method. In section [TV] we discuss how to solve numerically the governing transport equations by means 
of a lattice technique which bypasses the problem of solving separately the density and momentum equations. After 
validating the numerical method using some known results, we have made some applications to show its capability 
to predict the conductance and the volumetric flows in series of non uniform channels. In section |VT| we make some 
conclusive remarks and considerations. 



II. MODEL SYSTEM 



The model considered here treats the electrolyte solution as a ternary mixture. Positively and negatively charged 
hard spheres with point like charges located at their centers represent counter- and co-ions, while neutral spheres 
represents the solvent [33]. Each species, denoted by a, has mass m Q , diameter a aa , and charge z a e, e being the 
proton charge and z a the ion valence. The index a — 0, identifies the solvent, whose valence is zero, while a = ± 
identifies the two oppositely charged ionic species. The particles are immersed in a medium of constant and uniform 
dielectric permeability, e, and mutually interact via the potential 

U a Hr) = \ - < ~ aP (!) 
oo ll r < a p 



The non-coulombic part of the interaction is represented by an infinite repulsion which prevents the centers of any 
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two particles to get closer than the distance cr"' 9 = (cr aa + er ,3/3 )/2. We shall not consider other types of interactions 
such as attractive solvation forces between ions and solvent, which can be important to stabilize the solution |13| . 

In the kinetic description employed here, repulsive forces are accounted for by means of instantaneous collisions 
and treated a la Enskog, while the remaining forces are treated within the random phase approximation, disregarding 
correlations. Our treatment neglects the excess free energy stemming from electrostatic correlations, an approximation 
which can be justified assuming that it is small compared to the free energy excess associated with the excluded volume 
term. 

We first define the basic variables of our theory. Given the partial number densities n a (r,t), the total number 
density is 

n(r,i) = ]Tn Q (r,i), (2) 



the charge density p e (r 1 1) of the fluid is 



and the total mass density is 



p e (r,i) = e]>> Q n Q M), 

Oi 

Pm (r,t) = ^m Q n Q (r,i). 



(3) 
(4) 

(5) 



We also define the average velocities of each species, u a (r,t), and the average barycentric velocity as 



In the following, we neglect spatio-temporal variations of the temperature, T for the sake of simplicity. 



A. Continuum equations 

Before proceeding with the microscopic treatment, it is worth reviewing the macroscopic description of transport 
of charged liquids. Even the simplest theoretical description of electrokinetic phenomena requires the use of three 
physical models: a) the Poisson equation relating the electric field to its source, the spatial charge distribution; b) 
the constitutive Nernst-Planck equation giving the current of each ionic species as a function of fluid velocity, charge 
distribution and applied electric field; c) the hydrodynamic continuity and Navier-Stokes (NS) equations relating the 
density and fluid velocity to the body and pressure forces acting on the system. 

The three equations are coupled and their solution can be obtained by standard numerical methods [§]. Due to 
its low computational cost, the Poisson-Nernst-Planck equation represents the workhorse for diffusive electrokinetic 
phenomena |7J. It treats the ions as point particles mutually interacting via the Coulomb potential and experiencing 
a drag force stemming from the solvent. The solvent does not enter explicitly the description, but only determines 
the values of the dielectric constant and those of the transport coefficients. 

The ionic currents, J* = n ± u ± , have the following phenomenological Nernst-Planck expression 

3 ± (r,t) = -D ± Vn ± (r,t)-A ± ez ± n ± (r,i)V^(r,<) + n ± (r,t)u(r,t) (6) 

which is the superposition of three contributions: i) the one proportional to the gradient of the ionic concentration 
via the diffusion coefficient, D ± , ii) the migration current proportional to the gradient of the electrostatic potential 
ip(r,t) via the Nernst-Einstein mobility X ± = D ± /(ksT), on account of the friction exerted by the solvent on the 
ions, iii) and the convective current due to the displacement of the fluid with velocity u. In the Nernst-Planck 
approximation the cross coefficients that couple the diffusive transport of the different components are neglected. In 
the Stefan-Boltzmann approach, instead, such a coupling is preserved, but the determination of the cross-coefficients 
requires either a detailed microscopic theory or careful experimentation^. 

The electrostatic potential ip(r,t) is generated by the charge distribution p e (r,t) and by the fixed charges located 
on the pore surfaces and on the electrodes and satisfies the Poisson equation: 



vV(m) = -^M 



(7) 
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with boundary conditions — V^(r, t) = E(r)/e at the confining surfaces, where E(r) is the surface charge density. The 
Poisson-Nernst-Planck (PNP) equation is obtained by considering the local conservation law: 

^n ± (r,t)+V- J±(r,t) = (8) 

with J ± given by [6] and ip by the Poisson equation [7] To complete the picture one determines the velocity by 
considering the Navier-Stokes evolution equation for u(r,t): 

d t u + u • Vu = - — VP + + ^-V 2 u + I^±^V(V ■ u) (9) 

Pm Pm Pm Pm 

where P is the hydrostatic pressure, r\ and j]^ are the shear and bulk dynamic viscosities, respectively, and E the 
electric field which is the sum of the field due to the ions and of that due to the external charges. 

When considering a laminar flow in a straight slit-like channel, normal to the z direction, by imposing the vanishing 
of the current in the ^-direction intone obtains the Boltzmann ionic distribution along z: n ± (z) = c ± e -ez ± ip(z)/k B T ^ 
where ip( z ) is determined with the help of|7] The constants are fixed by boundary conditions, as discussed below. 

The stationary solutions of [9] of great interest are when the electrolytes are driven through narrow channels by an 
external electric field. Since the pioneering work of Smoluchowski on electro-osmotic phenomena, it is known that 
when the A^> << H the velocity is independent of the conduit size H. unlike pressure driven flows. This fact is very 
important and makes convenient to use electric driving to pump fluids in ultrafine capillaries. In nanometric channels 
this condition is hardly realized. 

For vanishing pressure gradients and under creeping flow conditions ( Stokes approximation ), one can use again 
equation [7] to eliminate the charge density p e , and obtain from [9] the following relation between the fluid velocity and 
the externally imposed electric field along the x direction: 

d 2 ifj(z) d 2 u x {z) 
^^^^^^ (10) 

with E x = — V x ^. By imposing the condition of vanishing velocity at the surfaces of the channel, located at z = 
and z — H, one finds the celebrated Smoluchowski relation: 

eK 

ux(*) = W*)-0— £ (ii) 
v 

where £ is the value of the potential at the shear plane where the velocity vanishes [5J. 

The main drawback of the macroscopic approach is that it holds only if the ions are treated as point-like particles 
migrating under the action of the self-induced electric field, diffusing in the solvent and convected by the flowing 
solvent. When transport occurs through narrow channels the deviation from this ideal gas representation becomes 
crucial. 



III. THE MICROSCOPIC APPROACH 



In the following, we shall introduce a microscopic representation of the electrolytic solution by switching to a kinetic 
description. We will extend the self-consistent dynamical method, that we recently employed to study neutral binary 
mixtures, to ternary charged mixtures. We consider the following set of Enskog-like equations governing the evolution 
of the one-particle phase space distributions f a (r,v,t) of each species: 

|f(r,v,f)+v Vf(r,v,t) + • ^/ a (r, v,t) = £ n Q "(r, v,t). (12) 

The left hand side of [12] represents the streaming contribution to the evolution, with F" being the external force 
acting on species a, and merely reflects the Liouville-like single particle dynamics, while the right hand side involves 
the interactions among the particles and will be the object of a series of approximations in order to render the theory 
numerically tractable. The first approximation |3 1] consists in separating the interaction into an harshly repulsive 
term, whose effect is represented as collisions between hard spheres, and terms whose gradients have a slower spatial 
variation and can be treated within a mean-field like approximation. According to the Enskog approach the effect of 
repulsion is obtained by considering pairs of particles dynamically uncorrelated before collisions, but with probabilities 
modulated by the static pair distribution function. The velocities of these pairs, however, emerge correlated after 



5 



a collision. The collision process is described within the revised Enskog theory (RET) proposed by van Beijeren 
and Ernst |22| which is known to give good results for the transport properties of hard spheres. Since even the 
RET operator f2 Q/3 (see Ref.[32 for an explicit representation) is awkward to manage, simplifying methods are very 
convenient: the RET operator is further simplified by invoking the approximation put forward by Dufty and coworkers 
|26U27| which separates the hydrodynamic from the non-hydrodynamic contributions to il a ^ and employs a projection 
technique. Extending such a prescription to the present case we obtain: 

^^(r,v,i)« 

P 

- u[r(r,v,t) - 0£(r, v,i)] + ■ (v - u(r,t))^(r,v,t) - ^V^(r) • ^/ a (r,v,t) (13) 

The term <& Q represents the sum of the internal forces of non electrostatic nature acting on species a. In addition, 
the electric potential ip is given by the solution of [7J Unlike the full collision operator of Enskog type, the last two 
terms in the right hand side of Eq. [13] act only on the hydrodynamic components of the distribution function f a . 
The non-hydrodynamic modes instead are driven by the first term in the right hand side of [13] which is a Bhatnagar- 
Gross-Krook (BGK) type of relaxation kernel. Its role is to drive the system towards local equilibrium over a time 
scale ui~ l , much shorter than the hydrodynamic time scales [34 , while the global equilibrium is controlled by the 
collisional force mentioned above. 

The BGK term contains the distributions functions and (j> a which have the following representations (see ref. 
|35j for details): 

4T (r, v, t) = n«(r, t)[^-]^ expf- ^^*^ (14) 
v V ' ' ' V ' n 2irk B T l y \ 2k B T J 1 ; 

and 

«(r,v,,) - f (r.v,«){. + '"°'"° (r -" - ■ <» - (15) 

Notice that in the case of a one-component fluid there is no difference between <frj_ and <fi a , since the velocities u a and 
u coincide and the standard BGK approximation involves the difference between the distribution f a and the local 
Maxwellian <j) a . The above prescription fullfils the indifferentiability principle which states that when all physical 
properties of the species are identical, the total distribution / = f a must obey the single species transport 
equation. 

The reason to use the modified distributions [15] instead ofPulis to obtain the correct mutual diffusion and hydrody- 



namic properties starting from 12 A simple BGK recipe, that is, by setting <f>^_ = 4> a , would lead to double counting 
the interactions on the diffusive properties. 

Equations ?? can be solved after specifying the forces F Q and the Coulombic force — eV-0(r) and the fluid 
velocities u and u a . 

From the knowledge of the / a 's it is possible to determine not only all the hydrodynamic fields of interest, but 
also the structure of the fluid at molecular scale. All equations which constitute the building blocks of the classical 
electrokinetic approach can be derived from the set of equations ??. In fact, the PNP and NS equations can be 
straightforwardly recovered by taking the appropriate velocity moments of the kinetic Enskog-Boltzmann equation. 
From the distributions, we compute the partial densities n a (r,t), 

n a (r,t) = Jdvf a (r,t) (16) 

the average velocities u a (r, t) of the species a, 

n a (r,t)u a (r,t)= J dvvf a (r,t), (17) 
and the kinetic contribution of component a to the pressure tensor, 

7rg(r,t) = m a /'dv(w i -^)(^--^-)r(r,v,f). (18) 



The explicit form of the internal forces <fr Q (r, t) featuring in 13 is determined by taking the first moment of the collision 



integral with respect to the velocity. By algebraic manipulations one observes that & a (r,t) is due to two-particle 
correlations and can be separated into three contributions [31,, as 

* Q (r,t) = F a > mf {r,t) +F a ' dra9 (r,t) +F a > msc (r,t). (19) 



6 



The first term arises from the gradient of the hard sphere chemical potential excess over the ideal gas value: 

F«.™/(r,t) = -V/C c (r,t) (20) 

with 

V^ xc = k B Tj2K P ) 2 [ dssg a0 (r,r + <j^s,t)n p (r + a al3 s,t) (21) 

13 J 

where the integration is over surface of the unit sphere and g a/3 {r, r + a a ^s, t) is the inhomogeneous pair correlation 
functions for hard spheres at contact (|r — r'| = cr a ^). The second term accounts for the frictional force between the 
ions and the solvent and depends linearly on their relative velocities: 

F a '* a9 (r,i) = -^ 7 Q,3 (r,t)(u Q (r,i)-i/(r,i)) (22) 



via the inhomogeneous friction tensor: 



where fi a ^ is the reduced mass /i a/3 = ^l,"^ for the colliding pair. Notice the local character of the relation between 
the drag force and the fluid velocity. Finally the particles experience a viscous force which has the following Enskog 
form: 

F a ' ci,c (r,t) =Y J 2(° a0 ) 2 y fJ ' a ' 3kBT [ dss 5Q(3 (r,r + cr Q,9 s,i)n' 3 (r + a a ' 3 s,t)s- (u^r + a^s) -u"(r)). (24) 

Since there is no exact theory of the pair correlation function in spatially inhomogeneous systems, the contact value 
of g a P is computed using a two-component generalization of the Fischer and Methfessel prescription[37 : one approxi- 
mates the inhomogeneous g a P by the corresponding bulk value evaluated when the partial densities take on the values 
of the smeared densities, n a (r): 

g^(r a ,r p ; |r° - r"| = o^) = ffSf Ifc ({f„(R^)}), (25) 

where the functions 

£ n (R^) = |^n a (R^)(cr aQ )" 

a 

are linear combinations of the smeared densities 

i e 



j- Q ( R «/3) = — / drn a (r - R° 
V a J v * 



over spheres of volume V a = n(a°" x ) 3 /6 centered at the point R"^ = r + r , where r a and r* 9 are the centers of 

the two spheres. The explicit expression of g^ lk is provided by an extension of the Carnahan-Starling equation to 
mixtures 



%uiMn}) - 1 _ ^ + 2 aa0 (1 _ 6)2 + 2 ( ^ J (1 _ 6)3 ■ (26) 

We remark that the present method requires the knowledge of the pair correlation function at contact, in contrast 
with the DDFT approach which requires instead the knowledge of free energy or equivalently the Ornstein-Zernike 
direct correlation function. However, in the spirit the method here employed is equivalent to the so called Weighted 
Density Approximation of equilibrium DFT. The latter does not contain information necessary to compute the trans- 
port coefficients, that in the present method are obtained from g a ° via the Enskog ansatz for the collision integrals. 
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A. Derivation of the equations of electrokinetics from the microscopic approach 



It is useful to show that in the limit of slowly varying fields the transport equations 12 reproduce the macroscopic 
evolution discussed in section [II A| To this purpose, we integrate [l2]w.r.t. the velocity and obtain the conservation 
law for the particle number of each species: 

§-n a (r,t) = -V • (n«M)(u Q (r,i) _ u(r,i))-V • (n Q (r,t)u(r,t)) = -V • J Q (r,i) , (27) 

In order to recover the PNP equation we consider the the momentum balance equation for the species a, which is 
obtained after multiplying by v and integrating w.r.t. v: 



n 



dt 



(r, t)uf(r, t)] + V, (n a (r, t)uf(v, t)uf(r, t) - n a (r, t)(uf(r, t) - Ul {v, t))(u?(r, t) - Uj (v, t)) 



- V< vK ' + -^^n a (r,t) + jV - n a (r,t) n a (r, i)V,^(r, t) (28) 

Using [20] we define the gradient of the global chemical potential of the individual species, fi a , as the sum of the ideal 
gas part and the excess part: 

n a (r,t)V lt , a (v,t) = VX,M)<% " n a (T,t)F?> mf (T,t) . (29) 

By assuming the existence of a steady current, we drop the non-linear terms in the velocities in the l.h.s. of |28| and 
neglect the viscous force contribution. The following approximated force balance is obtained 

V/i^M) - P ± (r) + ez ± V?A(r,i) « F ± ' drag (r, t). (30) 

The r.h.s. of |30| represents the drag force exerted on the particles of type a = ±, in reason of their different drift 
velocities. In dilute solutions the charged components are expected to experience a large friction arising only from 
the solvent while a negligible friction from the oppositely charged species, so that we further approximate 

F±.^(r,t) = - 7 ± (r,t)(u ± (r,t)-u(r,t)) (31) 



where we set u ss u°. In 31 the friction can be evaluated in uniform bulk conditions to be 



7* « ^2-Kk B T jtf \ g^n\a^f (32) 
3 V ui + m u 

with n° being the bulk density of the solvent and <?o± the bulk ion-solvent pair correlation function evaluated at 
contact. Finally, using eqs. [27] and [31] we obtain an expression for the ionic currents in terms of the microscopic 
parameters: 



Jf(r,t) = — ^-n ± (r,t)V/i ± (r,t) -ez ± n ± (r,tW^(r,t) + n ± (r, t)u(r, t) (33) 



which has the same form as the phenomenological Planck-Nernst current [6] with the full chemical potential gradient 
^ replacing the ideal gas chemical potential gradient, fcsTVm(n Q ) and 7 ± = ksT/D^ 1 and A* = l/7 ± . The total 
electric charge density current is 

± 7 

where the zero frequency electric conductivity a e \ is given by the Drude-Lorentz-like formula: 

a e; =e 2 (^n+M)+ M! n -(r,t)) (35) 

showing that the conductivity is due to collisions with the solvent and decreases as the solvent becomes denser (7 ± 
is an increasing function of n°) while increases with the number of charge carriers. It should be noticed that, in 
the approximation of a local friction, the Drude-Lorentz conductivity is local in space and time. However, in the 
numerical solution this local approximation is not imposed. 
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To derive the total momentum equation, already introduced on a phenomenological basis in [9] we sum |28| over the 
three components: 

d t u j {r,t)+u i (r,t)W l u :j (r,t) + —W l TT { i f ) + — ^ ez ± n ± (r,t)V j tp(r,t) 

Prn Pm _j_ 

-— n a (r 1 t)(F J a (r) + F?' mf (r,t) + F?' msc (r,t))=0. (36) 

Prn a=0,± 

where 7rj;^(r, i) = 7iy (r, i) is the total kinetic pressure which can be approximated as: 

ViTrjf } ~ %VjP w - r) (K) (iv^Vjtii + Vfo) , (37) 
with = ksT^2 a n a and 77 W = ^ Q ?i Q . Using the result of ref. |31) . we can write 

V n a (r, t)F a ' msc ( ri t) ~ -77^ V 2 u - {\r,^ + ^ C) )V(V ■ u). (38) 

a 

The non-ideal contribution to the shear viscosity is evaluated in uniform bulk conditions to be 

f,(°) = ~ £ V^V(^)V^X, (39) 



15 a/S 



30 



can be cast in the Navier-Stokes form 



while the bulk viscosity is rfc = %ff c ^ . In conclusion, 

1 O 71 ~T1 ~\~ T)h 

d t Uj + UiViUj = ViPSij - —Vjip + — VjVjUj + — -VjVjUi (40) 

Pm Pm Pm Pm 

with 77 = v {K) + V (C) and V 4 P - V.-Pm + £ a n Q (r, tJV^fr, t). 

Having established the correspondence between the macroscopic equations and those derived from the kinetic 
approach, we can now outline the method of solution. 

B. Entropic contribution to the Poisson-Boltzmann equation in the non primitive model 

Some recent studies on the non primitive model of ionic solution lead to a modified Poisson-Boltzmann equation 
|40l I42|. In the same spirit we derive within the present approach a similar equation. We consider only mixtures 



of equal sizes. At equilibrium, that is, when all the currents vanish, according to 33 the ionic densities satisfy the 
equations: 

fc B TVlnn ± (r,i) + V/w(r,£) +ez ± \7*jj(r) = (41) 
while the solvent satisfies the equation: 

fe B TVmn°(r,i)+V/w(r,f) =0 (42) 
Notice that since for equi-sized molecules /jf xc depends only on the overall local packing fraction one can write: 

fcsT^Vlmy^r,*) - Vln(n°(r,t))+ez ± V^(r) =0 (43) 

and by integrating one gets an expression relating the ionic densities to the solvent density 

n ± (r, t) = C±n\T, t) exp(- (44) 

where C is an integration constant to be found by normalization. 

As a result the electrostatic potential is coupled to the solvent density via the Poisson equation: 

VV(r, t) = - V(r, t) £ z±C± exp(- eE± *y r,t) ) (45) 

In the case of the restricted primitive model (no explicit solvent) the density of the "ghost" solvent is reabsorbed in 
the definition of the integration constant. The numerical test of such a relation is given below in section \V\ 
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IV. NUMERICAL METHOD 



The numerical solution of the coupled equations for the distribution functions eqs. 12 is achieved in the 
framework of the Lattice Boltzmann method, as presented in refs. [3~H [35] . The numerical method is a substantial 
modification of the conventional method used in fluid dynamics applications to the presence of hard sphere collisions 
|28| . In a nutshell, the Lattice Boltzmann method is based on the following steps. One discretizes the position 
coordinate, r, by introducing a Cartesian mesh whose lattice points are separated by a distance a. The continuous 
velocity, v, is also discretized by restricting its values to Q possible "states", v — > {c p } with p = 1,N. The discrete 
velocities are chosen as to connect neighboring spatial mesh points. For the present three-dimensional study, we use a 
19-speed lattice, consisting of one speed c p = (particle at rest on a mesh node), six discrete velocities with |c p | = 1 
and pointing towards the first mesh neighbors, and 12 particles with |c p | = \[2 pointing towards the second mesh 
neighbors. The continuous phase space distribution functions f a are replaced by the array f a (r,v,t) — > /5(r, t) and 
the velocity moments are evaluated as |5S] 

n a (r,t) = Y,fp (46) 

p 

and 

n Q M)u"(r,t) = ]Tc p / p «(r,i) (47) 

p 

The modification of the conventional Lattice Boltzmann method to the case of hard sphere dynamics is based on 
computing the forces encoded by eqs. [2T] [24] and [3l] via numerical quadratures. The presence of Coulomb forces 
requires the solution of the Poisson equation for the electrostatic potential generated by the mobile and surface 
charges. Its determination is achieved by using a successive over-relaxation method |14j . The speed of convergence 
of the Poisson solver is greatly enhanced by employing a Gauss-Siedel checker-board scheme in conjunction with 
Chebychev acceleration |14j . Neumann boundary conditions on the gradient of the electrostatic potential are imposed 
at the wall surface 

ft- Wires = (48) 
e 

where n is the normal to the surface, S. We impose an external voltage by emulating the presence of electrodes in 
real systems. The method implies imposing the value of the external electric field E x as the boundary condition 

VV- = -E x (49) 

at locations x = + and x = L~ . 

The electrostatic boundary conditions are imposed by using a second-order finite difference scheme, compatible with 
the overall accuracy of the Lattice Boltzmann method. The system is further taken to be periodic along the x direction 
for the fluid populations of each component. This choice implies that all hydrodynamic fields (density, currents, etc.) 
are periodic along the x direction. In order to accommodate the presence of the electrodes, the hydrodynamic fields 
exhibit a finite jump at the x = 0,L,2L, ... locations. Preliminary calculations have shown that this jump does not 
have major impact on the stationary quantities, in particular as regarding those in the pore region. No slip boundary 
conditions on the fluid velocities of each species are imposed by means of the bounce-back rule [28 . Once all forces 



on the right hand side of 12 are computed, we evolve the distribution function by a second-order accurate trapezoidal 
integration [36 . In addition, we compensate the harsh internal forces by auxiliary fields properly chosen. This method 
improves the accuracy and robustness of the numerical scheme without altering the structure and dynamics of the 
system. Full details can be found in ref. [36 



V. RESULTS 



We have conducted numerical experiments on channels of three different shapes, as illustrated in[T] Such choice is 
intended to capture the behavior of a large class of solid-state devices and biological ionic channels, i.e. macromolecular 
pores that control the passage of ions across cell membranes. As shown in[T] in system A the fluid mixture is confined 
to move between two uniformly charged parallel plates. Such a geometry is useful to study the electroosmotic flow 
and ionic current when the plate separation H becomes of the order of the Debye length [7] . In system B the surface 
charge of areal density £ covers only two rectangular regions of area A and longitudinal size a on opposite walls. 
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Finally, system C is obtained by considering a parallel slit containing a narrowing of length L. Such a striction forms 
an open wedge, where the planes are separated by a distance h at the small opening and H at the large opening. 
The wedge geometry recalls the shape of synthetic conical nanopores, which have been manufactured and studied in 
recent years |44fH5] in order to understand the origin of rectifying devices, as related to biological channels. 

It is important to remark that, for models B and C, away from the charged region we do not explicitly model a 
completely macroscopic reservoir. Rather, the far away region acts as a buffer that allows to recover electroneutrality 
sufficiently away from the charged region, as detailed in the following. Given the studied geometries, it is of little 
usage to compute characteristic curves and compare them directly with the experimental data obtained in electrolytic 
cells, as studied by Siwy and coworkers |441 145] . whereas our systems are more alike series of nano-devices. 

Throughout this paper we shall use lattice units defined in the following. Time is measured in units of the discrete 
step At = 1, lengths are measured in units of lattice spacing, the charge e and the mass m are assumed to be unitary. 
The thermal energy ksT is specified by fixing the thermal velocity whose value is vt = \fk^T~fm = l/\/3 in lattice 
units. The dielectric constant is obtained by fixing the Bjerrum length, Is = e 2 / (AireksT), which is the distance at 
which the thermal energy ksT is equal to the electrostatic energy between two unit charges. We remind that in water 
under ambient conditions and for two monovalent ions this length is Ib ~ 0.7 nm, therefore, we set Ib/& = 4. 

We consider mixtures of hard spheres of equi-sized diameter a aa = a = 4 in lattice units. The packing fraction is 
fixed to ^ J2 a n a (a aa ) 3 = 0.20. The surface charge density S is always taken to be negative and with typical values 
£cr 2 /e = 10 -3 . At initial time, the densities of the ions in the channel are set by fixing the Debye length in the bulk: 



Xd = J^f- = J^— (50) 

where n s is the average density of a single ionic species in the bulk. This relation is used to define the initial density 
of the minority ions (the coions) to n~ (r) = n s . Next, we impose a global electroneutrality condition in order to set 
the local density of counterions, as 



e / (n+(r) -n"(r))dr+ / E(r)dr = (51) 
Jv Js 

where the first integral extends over the volume of the system and the second over its surface. During the simulation, 
the number of particles of each species is fixed and, as time marches on, the densities rearrange according to their 
intrinsic evolutions. We fix the value at the Debye length to \d/& = 2.75, unless differently stated. 

To fix the physical units, by stipulating that a — 3 x 10~ 10 m, we consider a mixture with number densities 
n°a 3 — 0.392 and n ± cr 3 = 0.001408, corresponding to a 0.1 M solution. The friction coefficient 7 for such mixture is 
about 0.5 and the electric conductivity is er e / = 8 x 10~ 5 . Finally, the geometrical parameters of the pore are set to 
H/a = 12.5, h/a = 2.75 and a/a = 10. 

In the following we compare results at finite packing fraction with results obtained in the absence of HS collisions, 
that is, at negligible packing fraction. The former is useful to single out the role of molecular interactions while the 
latter corresponds to the simple Bhatnagar-Gross-Krook approximation, that is, the completely macroscopic picture. 

Although the numerical method produces complete time dependent solutions, we focus on studying stationary 
conditions for the sake of simplicity. For each simulated system, the convergence towards stationarity is fast (~ 2000 
iterations) and in general very stable. In the near future, we plan to explore time dependent phenomena, as in 
alternate electric field conditions. 

As a preliminary test, we validate the equilibrium properties of the system since in the absence of ionic currents, the 
theory predicts a simple local relation between the species densities and the electric potential obtained as immediate 
consequences of [44] and reading 

C_ 

n (r,tj C kb± 

and 

n + (r, t)n~(r, t) 



rag(r,t) 



= C+C- (53) 



that is, the ratio of the density profile probes the electrostatic field ip(r,t), whereas the product of the two profiles 
should be everywhere proportional to the square of the density of the solvent. These ratios are displayed inland 
[3] and the agreement is quite good for both tests. Joly and co-workers have performed a similar test in their study 
]40|. but in contrast to our method, they made the comparison by utilizing the density profiles obtained from MD 
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simulation of a ternary charged mixture. In the present calculation instead both the electric potential and the density 
profiles are obtained from the same numerical calculation in a self-consistent fashion. 

The usage of the full Poisson equation in connection with a hard sphere fluid disregards the effect of charge 
correlations [H]. This effect has been included by some authors [TT] using the Mean spherical approximation (MSA) 
represention of the direct pair correlation function. At the MSA level, the exclusion region of hard spheres is removed 
from the source term of [7J We have thus considered this correction to the electrostatic forces by removing the 
contribution 

lf = f dr'(p e (r')-Pe(r))^S (54) 



e J\r-r>\<a |r-r| J <^J\r-r'\<a I 1 " ~ l " |~ 

from the interaction. Notice that the second part of [54] is a simple rewriting of the left hand side that regularizes 
the integral and avoids numerical overflows. The profiles of the counterions with and without the correction [54] are 
reported in [4] It is apparent that at the molarities of interest as used in the present work, the effect of interactions 
arising from inner cores is very small for large and narrow pores and can be safely ignored in the simulations. 

A further test of the quality of the theoretical framework comes from comparing the equilibrium density profiles with 
the results of Molecular Dynamics simulations of a ternary mixture confined between two parallel slabs of opposite 
surface charges, that is, a slight variation of Model A. In the MD simulations the excluded volume interactions are 
modelled via a Lennard- Jones potential 



V LJ (r ij )=4e 



/ \ 12 




te) - 









(55) 



truncated at the cutoff, so that Vlj — for r%j > 2 1 / 6 er. The particles are repelled from the slab walls by the same 
type of interaction. In addition, the electrostatic forces are computed via the Ewald summation method for a periodic 
system (43]. The MD simulations are performed at constant volume and temperature, with parameters of e/fc^T = 2.5 
and a = 2. The comparison between the electrokinetic and MD results is done by determining the exclusion radius 
of the soft spheres used in MD, that is, by finding the effective distance r* such that Vlj(t*) — ksT. An analogous 
procedure is applied to determine the effective volume of the confining slab and, in this way, to match the densities of 
each component between the electrokinetic and MD simulations. The MD profiles are accumulated over 10 7 timesteps 
and for two different slabs of width H/cr — 10 and 4.5, filled with ~ 2000 and ~ 4000 particles resulting in packing 
fractions of 0.10 and 0.20, respectively. In order to accumulate sufficient statistics we perform the MD simulations 
at sufficiently high molarities corresponding to densities of n* 1 = 2.0 x 10~ 5 . The electrokinetic and MD profiles 
reported in [5] are in excellent agreement, in particular for the larger slab width. The quality of results does not seem 
to deteriorate with the packing fraction considered here. The above results are particularly reassuring in view of the 
small molarity accessible by the electrokinetic methodology, and overall provides good confidence on the simulation 
framework. 

One the of reasons for studying wedged channels having a nanometric orifice is the current rectification properties of 
this type of devices. In the past, such behavior has been ascribed to a ratchet mechanism arising from the joint effect 
of the asymmetric electrostatic potential within the pore region and departures from thermal equilibrium |44l I45j. 
Alternatively, it has been proposed that rectification arises from the inhomogeneous conductivity near the orifice that 
would produce asymmetry of fluxes of counterions and coions |48j . In this context, a third cause of rectification could 
be due to local charge accumulation in the wedge pore region, that is, an imperfect Donnan compensation of surface 
charges. It is thus interesting to look at the highly non-trivial distribution of neutral and charged species, as displayed 
in[6]for the contour plots for the density of species in model C for a surface charge of £ = —0.002. The plots show that 
the several walls present in the wedge geometry, strongly modulate the density profiles of each species as the effect 
of layering, and this layering overlaps with the double layer organization. We further notice that the concave corners 
act as strong accumulation points for the species, similarly to what shown some time ago by Dietrich and co-workers 
|47j . The reason for the accumulations is completely entropic and is such that the concave corners effectively attract 
particles, while convex corners exert repulsive forces. [6] reports the electric potential in the same geometry at zero 
and finite packing fraction and, surprisingly, the self-consistent solution for the potential does not depend sensibly on 
the presence of hard sphere interactions. 

We next evaluate the charge count along the x direction and in the absence of external voltage for model C. The 
charge count is defined as s ± (x) = j n (x, y, z)dydz, where S being the sectional area, which can be assimilated to 
a local surface charge of each species. The local excess charge As = (s + — s~) should partially compensate for the 
surface charge in the wedge region. [7] shows several interesting features that strongly depend on packing fraction. At 
first, we notice the appearance of oscillations in the profiles induced by the packing effects which are more pronounced 
near the narrow exit of the wedge. More importantly, the plots show that the compensation is imperfect, within the 
wedge region As stays around — £ and displays strong similarities between the zero and finite packing fraction cases. 
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[7] reveals how the molecular system departs from the assumption of local charge neutrality often utilized in electro- 
chemical studies. Only away from the pore region the excess charge decays to zero, that is, it restores electroneutrality. 
The longitudinal decay length is given by A^ on the side of the larger wedge mouth and shows a slightly more complex 
behavior on the other side. It is an interesting fact that the Debye length controls the longitudinal restoration of 
electroneutrality, perhaps not completely unexpected. In practical calculations, it is however important to determine 
the finite size effects, such as in the periodic system studied here. |8]shows that already channels of size L/\r> > 8 are 
sufficiently long to show a good convergence of the charge distribution. For our simulations, we have chosen L/a = 40 
as the reference channel extensions to obtain satisfactorily converged results. 

Once the voltage acts on the system, the charge unbalance modulates the charge organization in highly non-trivial 
ways. The charge distribution along the channel centerline of the wedged geometry is shown in|9]for different voltages. 
Counterions tend to accumulate and coions deplete in the wedge region for positive voltages. The situation reverses 
by reversing the bias and, under these conditions, we observe a strong charge piling up on the small mouth side. 
In [10 we turn our attention to the mass transport process and display the barycentric velocity profile, u z (x), of 



the fluid in the x direction in the case of model A. For the sake of comparison we also show the prediction of the 
Smoluchowski analytical formula 11 obtained by calculating the electrostatic potential ip(z) within the Debye-Hiickel 
linear approximation. We observe that the agreement between this formula and the results relative to the BGK system 
is very good, whereas the system with excluded volume interactions displays smaller velocities near the confining walls 
due to the larger value of the viscosity. 

Charge transport in the system is mostly due to ionic conduction rather than due to electro-osmotic transport. 
We therefore compute the number of mobile charge carriers in the pore region. As expected for model B , [TT] shows 
symmetric curves upon voltage inversion and rather flat curves. On the other hand, for the wedge system, the curves 
lose their symmetry and become more pronounced with the voltage. Only at the largest value of the surface charge 
density, one observes an appreciable deviation between the two channel shapes. In the lower panel of |11| the effect 
of the wedged geometry determines a non flat dependence of the number of the mobile charge carriers on the applied 
voltage. Also we notice that this number is lower in the presence rather than in the absence of HS collisions, since 
the excluded volume interaction hinders the pile up of charges in the wedged region. In any case, the undercharging 
of the pore at non-zero voltage is rather limited (< 20%). By the same token, the imperfect neutralization of the 
wedged region decreases the Drude-Lorentz conductivity for both positive and negative voltages, as compared to the 
zero voltage condition. The observed modest asymmetry in local charging also reflects in differences in conductivity 
for forward vs backward bias. However, rectification usually shows up as a much stronger effect |44l 145] . and therefore 
we exclude the fact that it is due to modulations of the occupation of the wedged region by mobile charge carriers. 

To demonstrate the bulk-like origin of electrical transport, in [12] we report the numerical data for the electric 
conductance, G, for the geometries of models B and C. The simulations provide the ionic current / and the potential 
difference V applied at the ends of the channel, so we computed the conductance according to the formula: 

I = GV. (56) 

An estimate of G can be derived using the expression G = , where S is the effective section of the channel, and 
given the geometry under consideration, S/L = 200/80 = 2.5. 

[13] reports the streaming current, that is, the ionic current arising from a pressure gradient Vp obtained for model 
B. The streaming current increases when increasing the Debye length, a behaviour that is expected by considering 
that the ionic current is proportional to the surface potential tp s , I — £ ^ 3 ^ Vp [3], and that the surface potential is 
roughly proportional to the Debye length. The presence of hard sphere collisions generates streaming currents smaller 
than without collisions, due to the increased value of the dynamic viscosity. 

It should be kept in mind that, in order to compare the computed conductances arising from either a electric 
field or a pressure gradient with the experimentally measured ones, the simulated system is the result of a periodic 
series of charged and uncharged regions of the pore. Nevertheless, the simulation data are important as they show the 
intimate nature of conduction. As a matter of fact, as the bulk concentration decreases, the excess surface charges due 
to the presence of the negatively charged surface of the channel becomes t he le ading contribution to the conduction 
mechanism. For very small values of n s the bulk formula a e i — 2^-n s (see 35 1 does not hold and must be modified 
in order to account for the contribution due to the excess charge arising from the counterions within the pore. By 



making the hypothesis that electroneutrality is not a global constraint (see 51 1 but rather holds locally, we simply 



replace n s — > n s — E/eiJ. The consequence of such replacement is the saturation of the conductance at low values 



of the electrolyte concentration n si as reproduced qualitatively in 12 Finally, the operational value of A^/cr = 2.75 



presently studied corresponds to being relatively far from the saturation regime. 
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VI. CONCLUSIONS 

We have studied nanometric electro-osmotic flows by means of a self-consistent kinetic description of the phase 
space distribution functions. The ternary mixture is composed by two oppositely charged species and a third neutral 
species representing the solvent, and is subject to external forces under non uniform conditions. The short-range 
interaction between molecules has been modeled by means of a hard sphere repulsion term and treated within the 
Revised Enskog theory. The Coulomb interaction has been treated within the mean field approximation, that is, by 
disregarding charge induced correlations. 

The presented theory was shown to be consistent with the phenomenological equations of macroscopic electrokinet- 
ics, such as the Smoluchowski and Planck-Nernst-Poisson equations. In essence, the approach bridges the microscopic 
with the hydrodynamic level and allows to incorporate the inhomogeneity of the ionic solution determined by the 
presence of the confining surfaces. 

In order to solve the system of equations for the multicomponent system, we have adopted a particular version of 
the Lattice Boltzmann method. We validated the theory and the numerical algorithm by considering the mass and 
charge transport under the effect of an applied electric fields in different channels having non uniform shapes and non 
uniform charging on the wall. 

In view of understanding the phenomenon of current rectification, as observed in wedged channels having an orifice 
of nanometric size, we have studied the departures from electroneutrality at molecular scale. We have found that 
electroneutrality is restored along the longitudinal direction within a Debye length. In addition, we observed mild 
undercharging of the wedged region at finite applied voltage. The asymmetry of the curves for positive vs negative 
voltage is modest, ruling out the possibility that this is at the origin of current rectification. 

Besides the geometries presented here, the proposed methodology can be applied to a larger class of problems in 
electrofluidics. For instance, it is possible to include semipermeable membranes, that is, channels impermeable to 
one type of ion and to the solvent. Other kinds of confining walls can be analyzed, such as chemically heterogeneous 
surfaces. A more comprehensive investigation of the role of differences in mass and size of the particles is currently 
underway. 
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FIG. 1: Geometries of the model systems A, B and C, with the shaded regions representing the wall. Each system contains 
three species of hard spheres of equal diameters a and masses m, but different charges: 0,1,-1, respectively. In model 
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iV+-iV _ = -2EwL/e, where N + and N~ are the number of positive and negative mobile charges, respectively, E is the surface 
charge density, and w the transversal dimension of the channel and L its length. Model B mimics a parallel slit channel, whose 
walls are decorated with negative surface charge density E in a region of length a. Charge neutrality , 7V + — iV - = -2Swa/e, 
holds. In Model C a wedge shaped region of length a connects two parallel slits. The surface charge sits only on the inclined 
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FIG. 2: Validation of the equilibrium conditions 52 and 53 for models A (upper) and B (lower panel). The left panels report 
the transversal profiles of the scaled potential eip/ksT (solid curves) and A + ln(n° /n + ) (circles) and B — ln(n°/n~) (squares), 
with A and B constants. The right panels report the transversal profile of the product n + n 7(n ) 2 . All profiles are taken at 
x = L/2. 
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FIG. 3: Same asp] for model C. The left panels report the transversal profiles of the scaled potential eip/ksT (solid curves) 
and A + ln(n°/n" (diamond symbols), with A a constant. The right panels report the transversal profile of the product 
n + n~ /(n ) 2 . The curves are transversal profiles taken at the pore larger mouth (black line and symbols) and at midpoint 
of the channel (red line and symbols). 
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FIG. 4: Model A: density profile of the counterions obtained with the full electrostatic interactions |7| and subtracted by 
the contribution of 54 for packing fraction 0.20, £ = —0.0003 and for channel width H/a — 12.5 (black line and red symbols) 
and H/a = 4.5 (red line and red symbols). The profiles have been rescaled by their maximum value for the sake of readability 
and are virtually indistinguishable. 
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FIG. 5: Comparison of density profiles as obtained from the electrokinetic method and Molecular Dynamic simulations for an 
equivalent soft sphere system. Simulation details are reported in the text. All profiles have been rescaled by its maximum value 
for the sake of readability. 
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FIG. 6: Model C: contours of scalar fields for the case with £ = —0.002 and voltage 0.3fcsT- Top to bottom: density of 
the neutral species, density of counterions, density of coions and electrostatic potential. The lowest panel is the electrostatic 
potential for the equivalent system without hard sphere forces. 
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FIG. 7: Model C: transversally averaged density profiles of the two ionic species with (upper panel) and without (lower panel) 
HS collisions. The black curve represents the counterion profile s + (x), the red curve the coion profile s~ , and the blue curve 
their difference (the local charge) . The latter has to be compared with the surface charge density profile (blue stepwise dashed 
curve). Finally, as a guide to to the eye the exponential decay with a characteristic Debye length, A_d, is reported as a violet 
dashed curve. 
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FIG. 8: Model C: density profiles along the channel centerline of counterions (upper panel) and coions (lower panel) without 
(left panel) and with (right panel) hard sphere collisions. Profiles for different channel lengths are shown: L = 120 (black), 
L = 160 (red) and L — 200 (blue) in lattice units. Data are for Scr 2 /e = —0.0032 and zero voltage. The dashed lines indicate 
the pore region. 
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FIG. 9: Model C: density profiles along the channel centerline of counterions (upper panel) and coions (lower panel) in the 
presence of hard sphere collisions. Data for voltages KjksT = —60 (black), —30 (red), (green), 30 (blue) and 60 (violet), 
taken at £<r 2 /e = —0.016 and L/a = 40. The dashed lines indicate the limits of the pore region. 
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FIG. 10: Model A. Barycentric velocity divided by its maximum value. The profile refers to an electrolyte mixture under 
a voltage V/ksT — 30 in the absence (squares) and in the presence of HS collisions (circles). The Smoluchowski analytical 
solution is also displayed (dashed curve). 



25 



0.1 



0.05 



-60 



-30 





Voltage/k T 



30 



60 



0.05 




-*=-*! 



IIJJ 



-60 



-30 





Voltage/k R T 



30 



60 



FIG. 11: Model C. Number of charge carriers in the pore region N c h — J v d i r[n + (r) + n (r)] for model B (upper panel) 

and model C (lower panel). The curves correspond to Ecr 2 /e = —0.0016 (circles), —0.0048 (squares), and —0.016 (triangles), 
respectively. Solid lines correspond are in the presence of HS collisions while the dashed lines are without HS collisions 
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FIG. 12: Double logarithmic plot of the conductance, G, versus number density of carriers n a for model B (open symbols) and 
model C (filled symbols), in the presence (circles) and in the absence (squares) of HS collisions. The conductance is reported 
in units of Go = e 2 /vt<t. The inset reports the conductance versus the Debye length. 
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FIG. 13: Double logarithmic plot of the streaming current versus number density of carriers n s for model B in the presence 
(circles) and in the absence (squares) of HS collisions. The applied force is 8.33 x 1CP 7 in units of v^/a. The ionic current is 
reported in units of Iq = evr/cr. The inset reports the current versus the Debye length. 



